drawing
Matthew Abbott / The New York Times


Introdução

Os incêndios florestais na Austrália são uma ocorrência regular que teve um papel significativo na formação da natureza do continente ao longo de milhões de anos. No entanto, os incêndios podem causar danos materiais significativos e levar à perda de vidas humanas e animais. Os incêndios florestais mataram cerca de 800 pessoas e bilhões de animais desde 1851. Estima-se que entre 2019 e 2020 os incêndios mataram pelo menos 33 pessoas e mais de 3 bilhões de animais.

Geralmente, os incêndios mais devastadores são precedidos por altas temperaturas, baixa umidade relativa e fortes ventos, que criam as condições ideais para a rápida propagação do fogo. A chuva é a maior aliada ao combate ao fogo e na mitigação do surgimento de novos de focos de incêndios. Outro problema que agrava a situação na Austrália, é a forte dependência do país de civis para conter seus incêndios, cerca de 90% dos bombeiros combatendo fogos florestais são voluntários.

Portanto, considerando os recursos limitados para combate aos incêndios florestais, ainda mais com o aumento das mudanças climáticas que estão tornando esses eventos mais frequentes e devastadores, saber se irá chover em uma região específica dadas algumas medições é importante para direcionar os esforços onde eles são mais necessários.

Dataset

O conjunto de dados utilizado contém cerca de 10 anos de observações diárias tiradas de estações meteorológicas em 49 locais na Austrália. O dataset é composto por informações como Data e Local da observação, dados meteorológicos do dia no local como Direção do vento, humidade e temperatura além de duas features categóricas: RainToday e RainTomorrow que indicam se choveu 1 mm ou mais no dia da observação e no dia seguinte, respectivamente.

Visualizando os Dados

dataset = tibble(read.csv('weatherAUS.csv'))
print(dataset)
## # A tibble: 145,460 x 23
##    Date       Location MinTemp MaxTemp Rainfall Evaporation Sunshine WindGustDir
##    <chr>      <chr>      <dbl>   <dbl>    <dbl>       <dbl>    <dbl> <chr>      
##  1 2008-12-01 Albury      13.4    22.9      0.6          NA       NA W          
##  2 2008-12-02 Albury       7.4    25.1      0            NA       NA WNW        
##  3 2008-12-03 Albury      12.9    25.7      0            NA       NA WSW        
##  4 2008-12-04 Albury       9.2    28        0            NA       NA NE         
##  5 2008-12-05 Albury      17.5    32.3      1            NA       NA W          
##  6 2008-12-06 Albury      14.6    29.7      0.2          NA       NA WNW        
##  7 2008-12-07 Albury      14.3    25        0            NA       NA W          
##  8 2008-12-08 Albury       7.7    26.7      0            NA       NA W          
##  9 2008-12-09 Albury       9.7    31.9      0            NA       NA NNW        
## 10 2008-12-10 Albury      13.1    30.1      1.4          NA       NA W          
## # ... with 145,450 more rows, and 15 more variables: WindGustSpeed <int>,
## #   WindDir9am <chr>, WindDir3pm <chr>, WindSpeed9am <int>, WindSpeed3pm <int>,
## #   Humidity9am <int>, Humidity3pm <int>, Pressure9am <dbl>, Pressure3pm <dbl>,
## #   Cloud9am <int>, Cloud3pm <int>, Temp9am <dbl>, Temp3pm <dbl>,
## #   RainToday <chr>, RainTomorrow <chr>

Usando a API do Google Maps para salvar as coordenadas geográficas de cada local de medição, como existem vários lugares com o mesmo nome em diferentes países anglófonos, é necessário concatenar o nome do país no final do local para garantir que as coordenadas estarão coerentes.

LocCount <- dataset %>% count(Location)

LocCount$Category <- paste(LocCount$Location, ", Australia")# Adding Australia to the end of location name since there are many cities with the same names in other countries

cities_df <- as.data.frame(LocCount)
locations_df <- mutate_geocode(cities_df, Category) # Getting the coordinates with the google maps API
locations <- as_tibble(locations_df)

print(locations)
## # A tibble: 49 x 5
##    Location          n Category                    lon   lat
##    <chr>         <int> <chr>                     <dbl> <dbl>
##  1 Adelaide       3193 Adelaide , Australia       139. -34.9
##  2 Albany         3040 Albany , Australia         118. -35.0
##  3 Albury         3040 Albury , Australia         147. -36.1
##  4 AliceSprings   3040 AliceSprings , Australia   134. -23.7
##  5 BadgerysCreek  3009 BadgerysCreek , Australia  151. -33.9
##  6 Ballarat       3040 Ballarat , Australia       144. -37.6
##  7 Bendigo        3040 Bendigo , Australia        144. -36.8
##  8 Brisbane       3193 Brisbane , Australia       153. -27.5
##  9 Cairns         3040 Cairns , Australia         146. -16.9
## 10 Canberra       3436 Canberra , Australia       149. -35.3
## # ... with 39 more rows
df <- merge(x = dataset, y = locations, by = "Location", all = TRUE) # Outer join to combine the orignial df with the location coordinates
# Creating Map Canvas
aus_coord <- geocode("Australia")

map <- get_googlemap(center = c(aus_coord$lon, aus_coord$lat), zoom = 4)

Mapa com as Médias de Temperaturas Máximas de Cada Local

MTemp <- df[!is.na(df$MaxTemp), ] %>% # Removing missing data
  select(Location, MaxTemp, lon, lat)


aggTemp <- aggregate(MTemp$MaxTemp, by = list(Locaton=MTemp$Location , Longitude=MTemp$lon , Latitude=MTemp$lat), FUN=mean) # Aggregating by Location and the Average  MaxTemp


ggmap(map) +
  geom_point(data = aggTemp,
             aes(x = Longitude, y = Latitude, size = x),
             color = "red", alpha = 0.3)  + 
  theme(legend.position = "left") +
  xlab("Longitude") + 
  ylab("Latitude") +
  labs(size="Average Max. Temperature (°C)")

Mapa de Calor de Pluvialidade

RF <- df[!is.na(df$Rainfall), ] %>% # Removing missing data
  select(Location, Rainfall, lon, lat)

aggRain <- aggregate(RF$Rainfall, by = list(Locaton=RF$Location , Longitude=RF$lon , Latitude=RF$lat), FUN=sum) # Aggregating by Location and the sum of Rainfall


YlOrBr <- c("#FFFFD4", "#FED98E", "#FE9929", "#D95F0E", "#993404")

ggmap(map) +
  stat_density_2d(
    data = aggRain,
    aes(x = Longitude, y = Latitude, fill = stat(x)),
    alpha = .1,
    bins = 30,
    geom = "polygon"
  ) +
  scale_fill_gradientn(colors = YlOrBr) +
  theme(legend.position = "left") +
  xlab("Longitude") + 
  ylab("Latitude") +
  labs(fill="Rainfall in mm")

Média de Precipitação por Mês

RD <- df[!is.na(df$Rainfall), ] %>%   # Removing missing data
  select(Date, Rainfall)

# substring(RD$Date,6,7) # Month only
RainMonth = aggregate(RD$Rainfall, by=list(Month=substring(RD$Date,6,7)), FUN=mean) # Average daily Rainfall by month


ggplot(data=RainMonth, aes(x=Month, y=x, group=1)) +
  geom_line()+
  geom_point() +
  xlab("Month") +
  ylab("Average Rainfall (mm per day)")

O gráfico acima demonstra que os meses de verão tem uma média de chuva muito maior que os meses de inverno, que tendem a ser mais secos.

Relação entre incidência solar e chuva

SunRain <- df[!is.na(df$Sunshine) & !is.na(df$Rainfall),] %>%   # Removing missing data
  select(Sunshine, Rainfall)


RainSun = aggregate(SunRain$Rainfall, by=list(Sunshine=SunRain$Sunshine), FUN=mean) # Average daily Rainfall by hours of sunshine


ggplot(data=RainSun, aes(x=Sunshine, y=x)) +
  geom_line()+
  xlab("Sunshine (hours in day)") +
  ylab("Average Rainfall (mm)")

Como é possível notar, existe uma relação inversa entre a incidência solar e a quantidade de chuva.

Incidência Solar por Mês

SunMonth <- df[!is.na(df$Sunshine) & !is.na(df$Date),] %>%   # Removing missing data
  select(Date, Sunshine)

MonthSun = aggregate(SunMonth$Sunshine, by=list(Month=substring(SunMonth$Date,6,7)), FUN=sum) # sunlight by month



ggplot(data=MonthSun, aes(x=Month, y=x, group=1)) +
  geom_line()+
  geom_point() +
  ylab("Sunshine (hours)") +
  xlab("Month")

Comparando o gráfico acima com o de Média de Precipitação por Mês, é reforçada a hipótese da relação inversa entre a incidência solar e a quantidade de chuva,

Precipitação em relação à direção do vento

DD<- df[!is.na(df$Rainfall), ] %>%  # Removing missing data
  select(Date, Rainfall,WindDir3pm, WindDir9am)

# Aggregating by Wind Direction and the avg Rainfall for 3pm and 9am
aggDir3 <- aggregate(DD$Rainfall, by = list(Dir3pm=DD$WindDir3pm), FUN=mean) 
aggDir9 <- aggregate(DD$Rainfall, by = list(Dir9am=DD$WindDir9am), FUN=mean)

# Transposing the Dataframes so the columns are the mean Rainfall for each direction and the rows are the differente times
final_df3 <- as.data.frame(t(aggDir3))
names(final_df3) <- lapply(final_df3[1, ], as.character)
final_df3 <- final_df3[-1,] 
rownames(final_df3) <- "3pm"

final_df9 <- as.data.frame(t(aggDir9))
names(final_df9) <- lapply(final_df9[1, ], as.character)
final_df9 <- final_df9[-1,] 
rownames(final_df9) <- "9am"

final_df <- rbind(final_df3, final_df9)

final_df <- final_df %>% mutate_if(is.character,as.numeric)

# Adding a Min and Max rows to the Dataframe to scale the chart 
max_min <- data.frame(
  N = c(4, 0), NNW = c(4, 0), NW = c(4, 0),
  WNW = c(4, 0), W = c(4, 0), WSW = c(4, 0),
  SW = c(4, 0), SSW = c(4, 0), S = c(4, 0),SSE = c(4, 0),
  SE = c(4, 0), ESE = c(4, 0),E = c(4, 0),
  ENE = c(4, 0), NE = c(4, 0),   NNE = c(4, 0)
)

rownames(max_min) <- c("Max", "Min")
final_df <- rbind(max_min, final_df)
op <- par(mar = c(1, 2, 2, 2))

radarchart(
    final_df, axistype = 1,
    color = c("#00AFBB","#FC4E07"), pfcol = scales::alpha(c("#00AFBB","#FC4E07"), 0.5),
    plwd = 2, plty = 1,cglcol = "grey", cglty = 1, cglwd = 0.8,    axislabcol = "grey", 
    caxislabels = c(0, 1, 2, 3, 4),vlcex = 0.7, vlabels = colnames(final_df),
    title = "Average Rainfall in mm by Wind Direction"
  )

legend(
  x = "bottomleft", legend = rownames(final_df[-c(1,2),]), horiz = TRUE,
  bty = "n", pch = 20 , col = c("#00AFBB", "#FC4E07"),
  text.col = "black", cex = , pt.cex = 1.5
  )

par(op)

Criando o modelo

Removendo colunas com muitos dados faltando e depois removendo entradas incompletas.

print(colMeans(is.na(dataset))) # Percentage of missing data in each column
##          Date      Location       MinTemp       MaxTemp      Rainfall 
##    0.00000000    0.00000000    0.01020899    0.00866905    0.02241853 
##   Evaporation      Sunshine   WindGustDir WindGustSpeed    WindDir9am 
##    0.43166506    0.48009762    0.07098859    0.07055548    0.07263853 
##    WindDir3pm  WindSpeed9am  WindSpeed3pm   Humidity9am   Humidity3pm 
##    0.02906641    0.01214767    0.02105046    0.01824557    0.03098446 
##   Pressure9am   Pressure3pm      Cloud9am      Cloud3pm       Temp9am 
##    0.10356799    0.10331363    0.38421559    0.40807095    0.01214767 
##       Temp3pm     RainToday  RainTomorrow 
##    0.02481094    0.02241853    0.02245978
datasetClean <- subset(dataset, select=-c(Evaporation, Sunshine, Cloud9am, Cloud3pm, Date)) # Removing features with high NA precentage
datasetClean <- datasetClean[complete.cases(datasetClean), ] # Removing NAs

One-Hot encoding as features categóricas.

WD3<-dummy(datasetClean$WindDir3pm)
WD9<-dummy(datasetClean$WindDir9am)
WGD<-dummy(datasetClean$WindGustDir)
Loc<-dummy(datasetClean$Location)

colnames(WD3) <- paste("WindDir3pm", colnames(WD3), sep = "_")
colnames(WD9) <- paste("WindDir9am", colnames(WD9), sep = "_")
colnames(WGD) <- paste("WindGustDir", colnames(WGD), sep = "_")
colnames(Loc) <- paste("Location", colnames(Loc), sep = "_")

datasetD <- cbind(datasetClean, WD3)
datasetD <- cbind(datasetD, WD9)
datasetD <- cbind(datasetD, WGD)
datasetD <- cbind(datasetD, Loc)


datasetD$RainToday <- factor(datasetD$RainToday,
             levels = c('No', 'Yes'),
             labels = c(0, 1))

datasetD$RainTomorrow <- factor(datasetD$RainTomorrow,
                   levels = c('No', 'Yes'),
                   labels = c(0, 1))

df<-tibble(datasetD)

Separando em datasets de treino e teste.

set.seed(1337)
split = sample.split(df$RainTomorrow, SplitRatio = 0.8)
training_set = subset(df, split == TRUE)
test_set = subset(df, split == FALSE)

Treinando o classificador Random Forest.

y <- training_set$RainTomorrow
X <- subset(training_set, select=-c(RainTomorrow))

classifier = randomForest(x = X,
                          y = y,
                          ntree = 30)

classifier
## 
## Call:
##  randomForest(x = X, y = y, ntree = 30) 
##                Type of random forest: classification
##                      Number of trees: 30
## No. of variables tried at each split: 10
## 
##         OOB estimate of  error rate: 15.53%
## Confusion matrix:
##       0    1 class.error
## 0 66370 3955  0.05623889
## 1 10077 9937  0.50349755

Calculando a acurácia do modelo.

# Accuracy Score

y_pred = predict(classifier, newdata = subset(test_set, select=-c(RainTomorrow)))

accuracy = sum(y_pred == test_set$RainTomorrow) / nrow(test_set)
print(accuracy)
## [1] 0.8572061

Considerando que previsões meteorológicas tem uma acurácia de cerca de 90% para previsões em um período de cinco dias(https://scijinks.gov/forecast-reliability/), a acurácia do modelo de aproximadamente 86% pode ser considerada razoável,

Matriz de confusão.

# Confusion Matrix
cm = table(test_set$RainTomorrow, y_pred)
print(cm)
##    y_pred
##         0     1
##   0 16774   807
##   1  2418  2586

Calculando o precision e recall.

recall = cm[2:2,2]/(cm[2:2,1]+cm[2:2,2])
precision = cm[2:2,2]/(cm[0:1,2]+cm[2:2,2])

cat("Precision:", precision, "\n")
## Precision: 0.7621574
cat("Recall:", recall)
## Recall: 0.5167866

Cerca de metade das vezes que o chove modelo predisse que não iria chover. Porém ele acerta a grande maioria das vezes quando prediz que vai ocorrer chuva no dia seguinte, o que é melhor do que se ele predisse que fosse chover e não chovesse. Portanto o modelo ainda pode ser melhorado ou ser utilizado outro tipo de modelo para fazer as predições.